PM 2.5 Concentrations

plotclr = brewer.pal(9, "PuBuGn")
class = classIntervals(Obs_value, 
                       9, 
                       style = "fixed", 
                       fixedBreaks = seq(min(Obs_value),
                       max(Obs_value), 
                       length = 9+1))
colcode <- findColours(class,plotclr)
plot(Longitude, Latitude, pch=20, col = colcode, main = "PM2.5 Measumrments on Select Sates")
US(xlim = range(Longitude),
   ylim = range(Latitude),
   add = T,
   lwd = 1,
   col = "gray")

map = GetMap.bbox(lonR = range(Longitude),
                  latR = range(Latitude), 
                  size = c(640,640), 
                  maptype = "hybrid")
PlotOnStaticMap(map)
convert_points = LatLon2XY.centered(map,
                                    Latitude,
                                    Longitude)
points(convert_points$newX,
       convert_points$newY, 
       col = colcode, 
       pch = 19)

x = Longitude
y = Latitude
model = lm(log(Obs_value) ~ x + y + State)
surf_est_surface = interp(x,
                          y,
                          model$fitted.values)
plotvar = surf_est_surface$z[!is.na(surf_est_surface$z)]
class = classIntervals(plotvar,
                       9,
                       style = "fixed",
                       fixedBreaks = seq(min(plotvar),
                       max(plotvar),
                       length = 9+1))
colcode = findColours(class, plotclr)
image.plot(surf_est_surface$x, 
           surf_est_surface$y, 
           surf_est_surface$z, 
           col = plotclr,
           breaks = seq(1.5,2.5, length = 10),
           zlim = c(1.5,2.5),
           xlab = "X", 
           ylab = "Y",
           main = "Estimated spatial trend")
US(xlim = range(Longitude),
   ylim = range(Latitude),
   add = T,
   lwd = 1,
   col = "red")

g = ggplot(data, aes(log(Obs_value))) + geom_density(aes(group = State, fill = State), alpha = 0.4) + ggtitle("Estimated Density of PM2.5 Concentration") + theme_economist()
ggplotly()
n = nrow(data)
dist_matrix = matrix(0 , nrow = n , ncol = n)
diff_matrix = matrix(0, nrow = n , ncol = n)
coordinates = matrix(cbind(Longitude , Latitude), nrow = n , ncol=2)
dist_matrix = rdist(coordinates, coordinates)
for(i in 1:n){
        for(j in 1:n){
                diff_matrix[i,j] <- (1/2)*(model$residuals[i]-model$residuals[j])^2 
        }
}
plot(as.numeric(dist_matrix), 
     as.numeric(diff_matrix), 
     type = "p", 
     pch = 20,
     col = "black",
     xlab = "Distance",
     ylab = "Squared differences",
     main = "Variogram Cloud Plot for PM2.5 Concentrations")

bins = (0.1)*((1:10)^1.8)
dt = data.frame(cbind(x, y, (model$residuals)))
empirical_variogram = variogram(model$residuals ~ 1 ,
                                locations = ~ x + y, 
                                dt,
                                boundaries = bins)
plot(empirical_variogram,
     col = "black",
     type = "p",
     pch = 20,
     main = "Empirical Semi-Variogram")

directional_empirical_variogram = variogram(model$residuals ~ 1,
                                            locations = ~ x + y, 
                                            dt,
                                            alpha = c(0,45,90,135),
                                            boundaries = bins)
plot(directional_empirical_variogram, main = "Directional Empirical Semi-Variogram")

exp_veriogram_model = fit.variogram(empirical_variogram,
                                    vgm(psill = 0.03, "Exp", 1, 0.01),
                                    fit.method = 2)
sph_veriogram_model = fit.variogram(empirical_variogram,
                                    vgm(psill = 0.03, "Sph", 1, 0.01),
                                    fit.method = 2)
gau_veriogram_model = fit.variogram(empirical_variogram,
                                    vgm(psill = 0.03, "Gau", 1, 0.01),
                                    fit.method = 2)
mat_veriogram_model = fit.variogram(empirical_variogram,
                                    vgm(psill = 0.03, "Mat", 1, 0.01, kappa = 1),
                                    fit.method = 2)
plot(empirical_variogram, exp_veriogram_model, main = "Exponential Model")

plot(empirical_variogram, sph_veriogram_model, main = "Spherical Model")

plot(empirical_variogram, gau_veriogram_model, main = "Gaussian Model")

plot(empirical_variogram, mat_veriogram_model, main = "Matern Model")

LS0tCnRpdGxlOiAiQklPU1RBVF8xIgpvdXRwdXQ6CiAgaHRtbF9ub3RlYm9vazogZGVmYXVsdAogIGh0bWxfZG9jdW1lbnQ6IGRlZmF1bHQKLS0tCgpQTSAyLjUgQ29uY2VudHJhdGlvbnMKCmBgYHtyLCBlY2hvPUYsIG1lc3NhZ2U9Riwgd2FybmluZz1GfQpsaWJyYXJ5KGdncGxvdDIpCmxpYnJhcnkobWFwcykKbGlicmFyeShSZ29vZ2xlTWFwcykKbGlicmFyeShzcCkKbGlicmFyeShyZ2RhbCkKbGlicmFyeShmaWVsZHMpCmxpYnJhcnkoUkNvbG9yQnJld2VyKQpsaWJyYXJ5KGNsYXNzSW50KQpsaWJyYXJ5KGdzdGF0KQpsaWJyYXJ5KG1hbmlwdWxhdGUpCmxpYnJhcnkoTUFTUykKbGlicmFyeShmaWVsZHMpCmxpYnJhcnkoYWtpbWEpCmxpYnJhcnkoZ2VvUikKbGlicmFyeShnZ3RoZW1lcykKbGlicmFyeShwbG90bHkpCgpzZXR3ZCgiL1VzZXJzL2JlbmRlbmlzc2hhZmZlci9Cb3ggU3luYy9VTSBGYWxsIDIwMTYvQklPU1RBVCA2OTYvRGF0YSIpCmRhdGEgPC0gcmVhZC50YWJsZSgiUGFydGljdWxhdGVfbWF0dGVyX2F1ZzIyX21pZHdlc3QudHh0IiwgaGVhZGVyID0gVCkKYXR0YWNoKGRhdGEpCmBgYAoKYGBge3IsIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9CgpwbG90Y2xyID0gYnJld2VyLnBhbCg5LCAiUHVCdUduIikKCmNsYXNzID0gY2xhc3NJbnRlcnZhbHMoT2JzX3ZhbHVlLCAKICAgICAgICAgICAgICAgICAgICAgICA5LCAKICAgICAgICAgICAgICAgICAgICAgICBzdHlsZSA9ICJmaXhlZCIsIAogICAgICAgICAgICAgICAgICAgICAgIGZpeGVkQnJlYWtzID0gc2VxKG1pbihPYnNfdmFsdWUpLAogICAgICAgICAgICAgICAgICAgICAgIG1heChPYnNfdmFsdWUpLCAKICAgICAgICAgICAgICAgICAgICAgICBsZW5ndGggPSA5KzEpKQoKY29sY29kZSA8LSBmaW5kQ29sb3VycyhjbGFzcyxwbG90Y2xyKQoKcGxvdChMb25naXR1ZGUsIExhdGl0dWRlLCBwY2g9MjAsIGNvbCA9IGNvbGNvZGUsIG1haW4gPSAiUE0yLjUgTWVhc3Vtcm1lbnRzIG9uIFNlbGVjdCBTYXRlcyIpClVTKHhsaW0gPSByYW5nZShMb25naXR1ZGUpLAogICB5bGltID0gcmFuZ2UoTGF0aXR1ZGUpLAogICBhZGQgPSBULAogICBsd2QgPSAxLAogICBjb2wgPSAiZ3JheSIpCgptYXAgPSBHZXRNYXAuYmJveChsb25SID0gcmFuZ2UoTG9uZ2l0dWRlKSwKICAgICAgICAgICAgICAgICAgbGF0UiA9IHJhbmdlKExhdGl0dWRlKSwgCiAgICAgICAgICAgICAgICAgIHNpemUgPSBjKDY0MCw2NDApLCAKICAgICAgICAgICAgICAgICAgbWFwdHlwZSA9ICJoeWJyaWQiKQoKUGxvdE9uU3RhdGljTWFwKG1hcCkKCmNvbnZlcnRfcG9pbnRzID0gTGF0TG9uMlhZLmNlbnRlcmVkKG1hcCwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgTGF0aXR1ZGUsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIExvbmdpdHVkZSkKcG9pbnRzKGNvbnZlcnRfcG9pbnRzJG5ld1gsCiAgICAgICBjb252ZXJ0X3BvaW50cyRuZXdZLCAKICAgICAgIGNvbCA9IGNvbGNvZGUsIAogICAgICAgcGNoID0gMTkpCgoKCmBgYAoKYGBge3IsIG1lc3NhZ2U9RkFMU0Usd2FybmluZz1GQUxTRX0KCnggPSBMb25naXR1ZGUKeSA9IExhdGl0dWRlCgptb2RlbCA9IGxtKGxvZyhPYnNfdmFsdWUpIH4geCArIHkgKyBTdGF0ZSkKCnN1cmZfZXN0X3N1cmZhY2UgPSBpbnRlcnAoeCwKICAgICAgICAgICAgICAgICAgICAgICAgICB5LAogICAgICAgICAgICAgICAgICAgICAgICAgIG1vZGVsJGZpdHRlZC52YWx1ZXMpCgpwbG90dmFyID0gc3VyZl9lc3Rfc3VyZmFjZSR6WyFpcy5uYShzdXJmX2VzdF9zdXJmYWNlJHopXQoKY2xhc3MgPSBjbGFzc0ludGVydmFscyhwbG90dmFyLAogICAgICAgICAgICAgICAgICAgICAgIDksCiAgICAgICAgICAgICAgICAgICAgICAgc3R5bGUgPSAiZml4ZWQiLAogICAgICAgICAgICAgICAgICAgICAgIGZpeGVkQnJlYWtzID0gc2VxKG1pbihwbG90dmFyKSwKICAgICAgICAgICAgICAgICAgICAgICBtYXgocGxvdHZhciksCiAgICAgICAgICAgICAgICAgICAgICAgbGVuZ3RoID0gOSsxKSkKCmNvbGNvZGUgPSBmaW5kQ29sb3VycyhjbGFzcywgcGxvdGNscikKCmltYWdlLnBsb3Qoc3VyZl9lc3Rfc3VyZmFjZSR4LCAKICAgICAgICAgICBzdXJmX2VzdF9zdXJmYWNlJHksIAogICAgICAgICAgIHN1cmZfZXN0X3N1cmZhY2UkeiwgCiAgICAgICAgICAgY29sID0gcGxvdGNsciwKICAgICAgICAgICBicmVha3MgPSBzZXEoMS41LDIuNSwgbGVuZ3RoID0gMTApLAogICAgICAgICAgIHpsaW0gPSBjKDEuNSwyLjUpLAogICAgICAgICAgIHhsYWIgPSAiWCIsIAogICAgICAgICAgIHlsYWIgPSAiWSIsCiAgICAgICAgICAgbWFpbiA9ICJFc3RpbWF0ZWQgc3BhdGlhbCB0cmVuZCIpCgpVUyh4bGltID0gcmFuZ2UoTG9uZ2l0dWRlKSwKICAgeWxpbSA9IHJhbmdlKExhdGl0dWRlKSwKICAgYWRkID0gVCwKICAgbHdkID0gMSwKICAgY29sID0gInJlZCIpCgpgYGAKCmBgYHtyLCBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQoKCmcgPSBnZ3Bsb3QoZGF0YSwgYWVzKGxvZyhPYnNfdmFsdWUpKSkgKyBnZW9tX2RlbnNpdHkoYWVzKGdyb3VwID0gU3RhdGUsIGZpbGwgPSBTdGF0ZSksIGFscGhhID0gMC40KSArIGdndGl0bGUoIkVzdGltYXRlZCBEZW5zaXR5IG9mIFBNMi41IENvbmNlbnRyYXRpb24iKSArIHRoZW1lX2Vjb25vbWlzdCgpCmdncGxvdGx5KCkKCmBgYAoKYGBge3IsIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9CgpuID0gbnJvdyhkYXRhKQpkaXN0X21hdHJpeCA9IG1hdHJpeCgwICwgbnJvdyA9IG4gLCBuY29sID0gbikKZGlmZl9tYXRyaXggPSBtYXRyaXgoMCwgbnJvdyA9IG4gLCBuY29sID0gbikKCmNvb3JkaW5hdGVzID0gbWF0cml4KGNiaW5kKExvbmdpdHVkZSAsIExhdGl0dWRlKSwgbnJvdyA9IG4gLCBuY29sPTIpCgpkaXN0X21hdHJpeCA9IHJkaXN0KGNvb3JkaW5hdGVzLCBjb29yZGluYXRlcykKCmZvcihpIGluIDE6bil7CiAgICAgICAgZm9yKGogaW4gMTpuKXsKICAgICAgICAgICAgICAgIGRpZmZfbWF0cml4W2ksal0gPC0gKDEvMikqKG1vZGVsJHJlc2lkdWFsc1tpXS1tb2RlbCRyZXNpZHVhbHNbal0pXjIJCiAgICAgICAgfQp9CgoKCnBsb3QoYXMubnVtZXJpYyhkaXN0X21hdHJpeCksIAogICAgIGFzLm51bWVyaWMoZGlmZl9tYXRyaXgpLCAKICAgICB0eXBlID0gInAiLCAKICAgICBwY2ggPSAyMCwKICAgICBjb2wgPSAiYmxhY2siLAogICAgIHhsYWIgPSAiRGlzdGFuY2UiLAogICAgIHlsYWIgPSAiU3F1YXJlZCBkaWZmZXJlbmNlcyIsCiAgICAgbWFpbiA9ICJWYXJpb2dyYW0gQ2xvdWQgUGxvdCBmb3IgUE0yLjUgQ29uY2VudHJhdGlvbnMiKQpgYGAKCgpgYGB7cn0KCmJpbnMgPSAoMC4xKSooKDE6MTApXjEuOCkKZHQgPSBkYXRhLmZyYW1lKGNiaW5kKHgsIHksIChtb2RlbCRyZXNpZHVhbHMpKSkKCmVtcGlyaWNhbF92YXJpb2dyYW0gPSB2YXJpb2dyYW0obW9kZWwkcmVzaWR1YWxzIH4gMSAsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbG9jYXRpb25zID0gfiB4ICsgeSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZHQsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgYm91bmRhcmllcyA9IGJpbnMpCgpwbG90KGVtcGlyaWNhbF92YXJpb2dyYW0sCiAgICAgY29sID0gImJsYWNrIiwKICAgICB0eXBlID0gInAiLAogICAgIHBjaCA9IDIwLAogICAgIG1haW4gPSAiRW1waXJpY2FsIFNlbWktVmFyaW9ncmFtIikKCmBgYAoKCmBgYHtyfQpkaXJlY3Rpb25hbF9lbXBpcmljYWxfdmFyaW9ncmFtID0gdmFyaW9ncmFtKG1vZGVsJHJlc2lkdWFscyB+IDEsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbG9jYXRpb25zID0gfiB4ICsgeSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZHQsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgYWxwaGEgPSBjKDAsNDUsOTAsMTM1KSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBib3VuZGFyaWVzID0gYmlucykKCnBsb3QoZGlyZWN0aW9uYWxfZW1waXJpY2FsX3ZhcmlvZ3JhbSwgbWFpbiA9ICJEaXJlY3Rpb25hbCBFbXBpcmljYWwgU2VtaS1WYXJpb2dyYW0iKQoKYGBgCgoKYGBge3J9CgpleHBfdmVyaW9ncmFtX21vZGVsID0gZml0LnZhcmlvZ3JhbShlbXBpcmljYWxfdmFyaW9ncmFtLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB2Z20ocHNpbGwgPSAwLjAzLCAiRXhwIiwgMSwgMC4wMSksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGZpdC5tZXRob2QgPSAyKQoKc3BoX3ZlcmlvZ3JhbV9tb2RlbCA9IGZpdC52YXJpb2dyYW0oZW1waXJpY2FsX3ZhcmlvZ3JhbSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgdmdtKHBzaWxsID0gMC4wMywgIlNwaCIsIDEsIDAuMDEpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBmaXQubWV0aG9kID0gMikKCmdhdV92ZXJpb2dyYW1fbW9kZWwgPSBmaXQudmFyaW9ncmFtKGVtcGlyaWNhbF92YXJpb2dyYW0sCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHZnbShwc2lsbCA9IDAuMDMsICJHYXUiLCAxLCAwLjAxKSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZml0Lm1ldGhvZCA9IDIpCgptYXRfdmVyaW9ncmFtX21vZGVsID0gZml0LnZhcmlvZ3JhbShlbXBpcmljYWxfdmFyaW9ncmFtLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB2Z20ocHNpbGwgPSAwLjAzLCAiTWF0IiwgMSwgMC4wMSwga2FwcGEgPSAxKSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZml0Lm1ldGhvZCA9IDIpCgoKcGxvdChlbXBpcmljYWxfdmFyaW9ncmFtLCBleHBfdmVyaW9ncmFtX21vZGVsLCBtYWluID0gIkV4cG9uZW50aWFsIE1vZGVsIikKcGxvdChlbXBpcmljYWxfdmFyaW9ncmFtLCBzcGhfdmVyaW9ncmFtX21vZGVsLCBtYWluID0gIlNwaGVyaWNhbCBNb2RlbCIpCnBsb3QoZW1waXJpY2FsX3ZhcmlvZ3JhbSwgZ2F1X3ZlcmlvZ3JhbV9tb2RlbCwgbWFpbiA9ICJHYXVzc2lhbiBNb2RlbCIpCnBsb3QoZW1waXJpY2FsX3ZhcmlvZ3JhbSwgbWF0X3ZlcmlvZ3JhbV9tb2RlbCwgbWFpbiA9ICJNYXRlcm4gTW9kZWwiKQoKYGBgCgoK